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Abstract 

In the first part of this series of two papers, we proposed a theoretical formalism that enables one 
to model and categorize heterogeneous materials (media) via two-point correlation functions 5*2 and 
introduced an efficient heterogeneous-medium (re) construction algorithm called the "lattice-point" 
algorithm. Here we discuss the algorithmic details of the lattice-point procedure and an algorithm 
modification using surface optimization to further speed up the (re) construction process. The 
importance of the error tolerance, which indicates to what accuracy the media are (re)constructed, 
is also emphasized and discussed. We apply the algorithm to generate three-dimensional digitized 
realizations of a Fontainebleau sandstone and a boron carbide/aluminum composite from the two- 
dimensional tomographic images of their slices through the materials. To ascertain whether the 
information contained in 52 is sufficient to capture the salient structural features, we compute 
the two-point cluster functions of the media, which are superior signatures of the microstructure 
because they incorporate topological connectedness information. We also study the reconstruction 
of a binary laser-speckle pattern in two dimensions, in which the algorithm fails to reproduce 
the pattern accurately. We conclude that in general reconstructions using 52 only work well for 
heterogeneous materials with single-scale structures. However, two-point information via 52 is 
not sufficient to accurately model multi-scale random media. Moreover, we construct realizations 
of hypothetical materials with desired structural characteristics obtained by manipulating their 
two-point correlation functions. 

PACS numbers: 05.20.-y, 61.43.-j 
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I. INTRODUCTION 



Random heterogeneous multiphase materials or media are ubiquitous. Examples include 
composites, porous media, biological materials as well as cosmological structures, and their 
macroscopic properties are of great interest In the first part of this series of two papers 
[2! (henceforth referred to as paper I), we proposed a theoretical formalism to model and 
categorize heterogeneous materials via two-point correlation functions 5*2 (r), which can be 
interpreted as the probabihty of finding two points separated by the the dsiplacement vector 
r in one of the phases In particular, we introduced the idea of the two-point correlation 
function space and its basis functions. In general, 5*2 of a medium can be expressed by 
a map p on the associated basis functions, which is composed of convex-combination and 
product operations. We also suggested a set of basis functions by examining certain known 
realizable analytical two-point correlation functions. Moreover, we introduced an efficient 
isotropy-preserving S'2-sampling algorithm, namely the lattice-point algorithm, but left the 
details of the algorithm for another paper. 

Here we will provide algorithmic details of the lattice-point methodology and consider 
several nontrivial applications to illustrate the practical utility of our theoretical formalism. 
In particular, we will apply the lattice-point algorithm to generate three-dimensional (3D) 
digitized realizations of a Fontainebleau sandstone and a boron carbide/aluminum compos- 
ite from two-dimensional (2D) tomographic images of their slices through the materials. 
To justify whether the reconstructions are successful, one should also measure other sta- 
;istical descriptors of the media Here we compute the two-point cluster function C2 
^ (see definition and discussions in Sec. IV), which contains nontrivial "connectedness" 
information of the phases of interest. By comparing C2 of the target and reconstructed 
media, we demonstrate that S2 is indeed sufficient to capture the salient structural features 
in these cases. We also study the reconstruction of a binary laser-speckle pattern in two 
dimensions, and find that the algorithm fails to reproduce the target pattern accurately. 
We conclude that in general reconstructions using 5*2 only work well for heterogeneous 
materials with single-scale structures (those having only one characteristic length scale). 
However, two-point information via 5*2 is not sufficient to accurately model multi-scale me- 
dia (those composed of structural elements associated with multiple characteristic length 
scales). Moreover, we will show how the proper convex combinations of basis functions en- 
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able one to obtain two-point correlation functions with desired properties and thus enable 
one to generate a variety of structures with controllable morphological features. Note that 
here we will mainly focus on modeling heterogeneous materials in three dimensions, which 
complements the two-dimensional examples that we considered in paper I. 

Yeong and Torquato formulated the (re) construction problem as an energy-minimization 
problem using simulated annealing 5|. This has become a very popular (re) construction 
technique [g], Q, [sl, I^, l^- In this method, a nonnegative objective function E, called the 
"energy," is defined as the sum of squared differences between the target and sampled cor- 
relation functions. In our case, the energy is given by 

^ = E[^2(r.)-^2(r,)f, (1) 

i 

where S2{r) and 5*2 (r) is the target and sampled two-point correlation function, respectively. 
An important issue in the (re) constructions based on the method of simulated annealing is 
the choice of the energy threshold Eth, i.e., the error tolerance of discrepancies between 
the correlation functions of generated medium and the imposed ones. When the energy of 
the (re) constructed medium is below Eth, the (re) construction process is terminated. The 
energy threshold is a key indicator of how accurately the medium is (re) constructed. In 
our previous work, Eth was always chosen to be a very small number, e.g., 10~^; however, 
no quantitative analysis on why such a value should be chosen was given. In this paper, 
we will discuss in detail the significance of the energy threshold for both the orthogonal 
S'2-sampling algorithm (only sampling 5*2 along convenient orthogonal directions) {3,5] and 
the lattice-point algorithm. 

In the lattice-point procedure, we consider the digitized medium (pixel system) to be 
a "lattice-gas" system, in which the black pixels behave like nonoverlapping "particles" 
moving from one lattice site to another. The two-point correlation function of the medium 
is obtained by binning all the distances between the black pixels and dividing the number 
in each bin by the total number of distances between all lattice sites within that bin. To 
generate a trial configuration, a randomly selected "particle" (black pixel) is given a random 
displacement subjected to the nonoverlapping constraint. Only the distances between the 
moved "particle" and all the others need to be recomputed in order to obtain 5*2 of the trial 
configuration. In this way, all directions in the digitized medium are effectively sampled; 
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moreover, the complexity of the algorithm is linear in Nb, i.e., the number of nonoverlapping 
particles. 

For heterogeneous materials containing well-defined inclusions, we will show that the 
computational speed of the (re) construction process can be increased by an algorithm mod- 
ification using surface optimization, which is essentially a biased pixel-selection procedure. 
This algorithm modification is based on the fact that in the later stages of (re) constructions, 
further refinements of the configurations are achieved by the moves of black pixels on the 
surfaces of formed pixel-clusters. Thus, only the "surface pixels" are selected and given 
random moves to generate trial configurations, i.e., the surfaces are optimized. The phys- 
ical analog of this process is solidification, in the later stage of which the re-arrangements 
of solutes only occur on the surfaces of nucleated particles. We will see in the following 
that one can obtain a much lower error (discrepancies between the correlation functions of 
the (re) constructed medium and the imposed ones) when surface optimization is properly 
applied. This improvement enables one to quantitatively study the non-uniqueness problem 



of the (re) constructions 
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The rest of the paper is organized as follows: In Sec. II, we describe the lattice-point 
algorithm and the algorithm modification using surface optimization in great detail. The 
choice of different pixel-lattices is also discussed. In Sec. Ill, we discuss the significance 
of the energy threshold for both the orthogonal S'2-sampling method and the lattice-point 
algorithm. In Sec. IV, we suggest a general form of the map p and apply the theoretical 
formalism to model a Fontainebleau sandstone, a boron carbide/aluminum composite and 
a binary laser-speckle pattern. We also show how to construct materials with structural 
properties of interest by manipulating the parameters in the basis functions and choosing 
proper combination coefficients. In Sec. V, we make concluding remarks. 



II. ALGORITHMS FOR GENERATING HETEROGENEOUS MATERIALS 

In paper I, we derived the exact algebraic equations for the (re) construction problem and 
showed that the equations have an infinite number of solutions and cannot be solved rigor- 
ously. In principle, the Yeong-Torquato scheme enables one to obtain one of the solutions 
efficiently. The most time-consuming steps of the scheme are the samplings of the two-point 
correlation function at every trial configuration. An efficient S'2-sampling method would 
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dramatically speed up the (re) construction process. Furthermore, an isotropy-preserving 
algorithm is required in (re) constructions whenever a radial two-point correlation function 
5*2 (r) (r = |r|) is employed for the case of statistically homogeneous and isotropic media. 

A. Lattice-Point Algorithm 

The lattice-point algorithm is designed to sample the digitized representation of a statisti- 
cally homogeneous and isotropic medium in all possible directions efficiently. For simplicity, 
we will illustrate the idea in 2D. Implementation of the algorithm in three dimensions is 
a straightforward extension. Instead of considering the digitized medium as a collection 
of black and white pixels, one can think of the medium as a lattice-gas system: the black 
pixels are the "gas molecules" and the white pixels are unoccupied lattice sites, as shown in 
Fig. [H The "gas molecules" are free to move from one lattice site to another, subject to the 
impenetrability condition, i.e., each lattice site can only be occupied by one "gas molecule". 
Thus, the black pixels behave like hard particles and the volume fraction of black phase is 
conserved during the evolution of the system. 




FIG. 1: Digitized medium as lattice-gas system: (a) White pixels as unoccupied lattice sites, (b) 
Black pixels as non-overlapping "gas molecules". 

For a statistically homogeneous and isotropic particle system (e.g., an equilibrium hard- 
sphere system or an equilibrium lattice-gas system in (i-dimensional Euclidean space K*^), 
the most basic statistical descriptor of the spatial correlations of the particles is the pair 
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correlation function g2{r) The quantity pg2{f)si{r)dr is proportional to the conditional 
probability of finding the center of a particle in a d-dimensional spherical shell of volume 
Si{r)dr, given that there is another particle at the origin. Here p = N/V is the number 
density of the system and Si(r) is the surface area of a d-dimensional sphere of radius r, 
which is given by 



d-l 



where r(x) is the Euler-Gamma function. Hence, for a finite system, integrating pg2{r) over 
the volume yields (A^ — 1), i.e., all the particles except the one at the origin. Alternatively, 
pg2{r)si{r)dr is the average number of particles at a radial distance between r and r + dr 
from a reference particle. Practically, g2{r) can be obtained from simulations by generating 
a histogram for the number of particles n{r) contained in a concentric shell of finite thickness 
("bin" width) Ar at radial distance r from a arbitrarily chosen reference particle. 

The two-point correlation function 5*2 (r) of a statistically homogeneous and isotropic 
medium can be interpreted as the probability that both of ends of a randomly oriented line 
segment with length r fall into the phase of interest, say, the "black" phase. By comparison, 
we can see that the method of computing g2{r) of an isotropic particle system implies a 
natural way of obtaining 5*2 (r) of a statistically isotropic digitized medium, which efficiently 
uses all possible vector-information in the medium. For each black pixel (or "gas molecule" ) 
i, the distances between pixel i and all the other pixels j are computed and binned to 
generate a histogram for the number of black pixels separated from each other by distance 
r. Another histogram for the number of lattice sites separated from a reference site by 
distance r is also generated by computing and binning all possible site-separation distances. 
Suppose the histograms for black pixels and lattice sites are stored in the array BN[r] 
and SN[r], respectively. It is easy to show that the two-point correlation function can be 
obtained from 



S2[r] = BN[r]/SN[r], (3) 

(see Fig. [2]). In other words, to compute S2 we first calculate the fraction of occupied lattice 
sites separated by distance r (r^ is an integer) from a reference black pixel. Then we average 
this fraction over all black pixels (by choosing every black pixel as reference pixel once) to 
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obtain 5*2 (r). This procedure is consistent with the geometrical probabihty interpretation 
of 5*2. For the digitized media, a natural bin width could be the characteristic size of the 
pixel of the lattice, e.g., the edge length of a square pixel if square lattice is used. 




(a) (b) 

FIG. 2: (color online). Centers of particles contained in the concentric shell of a reference particle: 
(a) In a equilibrium hard-disk system, (b) In a lattice-gas system (or a digitized medium). 

At each step in the simulated annealing process, a trial configuration is generated by 
moving a randomly selected black pixel to an unoccupied lattice site. A configuration array 
can be used to speed up this process. In our 2D implementation, the configuration array is 
a 2D array with the entries being 1 and 0, corresponding to the occupied and unoccupied 
lattice sites, respectively. When a trial move is made, first we need to test whether it violates 
the impenetrability condition by checking whether the underlying lattice site is occupied or 
not. This process just requires constant access time to an entry of the configuration array. 
Note that several trial moves can be attempted before a trial configuration is found for a 
high-density system. After a trial configuration is generated, the old position information 
of the selected pixel is stored in the array designated PO. 

The next step is to recompute the two-point correlation function for the trial configura- 
tion. To do this efficiently, a distance matrix is set up when the system is initialized and 
updated if a trial configuration is accepted. Note that this matrix is symmetric. For each 
trial configuration, only the position of a randomly selected black pixel is changed, e.g., the 
k^^ pixel. Thus, only the k^^ row and column of the distance matrix need to be updated, 
which requires Nb operations {Nb is the total number of black pixels in the system). The 
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old entries in the k^^ row (or column) of the distance matrix are stored in the array DO. 
It is also unnecessary to recompute the entire array BN[r]. Recall that BN[r] contains the 
binned number of black pixels separated from each other by distances between r and r + Ar 
(Ar is the bin width). Once the k^^ pixel is selected, the distances between pixel k and all 
the other pixels are binned and stored in the array BNO[r]. After the trial configuration is 
generated, the new distances between pixel k and all the other pixels are binned and stored 
in the array BNN[r]. The array BN[r] is updated as follows for each r: 



5*2 [r] is then recomputed using Eq. (j3]) and the trial configuration is accepted with the 
probability given by Eq. (57) in paper I. If the trial configuration is rejected, all information 
of the old configuration can be restored easily using the arrays PO, DO, BNO[r] and 
BNN[r]. For example, BN[r] of the old configuration can be restored by 



The above procedures are repeated until the energy of the (re) constructed medium is 
below the energy threshold [see Eq. ([T]) in Sec. Ill] or the total number of evolutions reaches 
the prescribed limit value. 

B. Algorithm Modification Using Surface Optimization 

In the (re) construction of media composed of well-defined "particles" or large clusters, we 
find out that in the later stages of simulated annealing, the random pixel-selection process 
is very inefficient. Many trial moves are rejected because a majority of selected pixels are 
inside the formed "particles" or clusters and it is energetically unfavorable to move them 
outside. This requires the use of a biased pixel-selection process, i.e., surface optimization. 

The idea of surface optimization is analogous to the physical process of solidification: 
when the nuclei of proper sizes have been formed, they capture more solutes from surrounding 
solution to further decrease the total free energy. In the simulated annealing process, once the 
"nuclei" are formed, the random pixel-selection process is replaced by a biased pixel-selection 
process, i.e., only those in the surrounding "solution" or on the surface of a "nucleus" are 
selected to be moved to or along the surface of randomly selected "nuclei" . 



BN[r] - BNO[r] + BNN[r] BN[r]. 



(4) 



BN[r] - BNN[r] + BNO[r] BN[r]. 



(5) 
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It is natural and easy to incorporate the surface optimization with the lattice-point al- 
gorithm because the black pixels are already considered as "molecules". Each black pixel 
is assigned a "free energy", which is the minus of the number of the nearest neighbors of 
that pixel. If a pixel is inside a "nuclei" , it has the largest number of nearest neighbors and 
the lowest "free energy", which equals —4 if a square lattice is used and —6 if a triangular 
lattice is used. If the pixel is on the surface of a "nucleus" or in the surrounding "solution" , 
its "free energy" will be relatively higher. The highest free energy is 0, which means the 
pixel is separated from all the others. Thus, the black pixels are grouped into two subsets: 
the low-energy subset that contains pixels with lowest "free energy" and the high-energy 
subset that contains the other pixels. When the trial move is made, only the pixels in the 
high-energy subset are selected to bias the move. If the trial move is accepted, the free 
energy of the moved pixel and its neighbor pixels are recomputed; the subsets are updated. 




50 100 150 200 

n 



FIG. 3: Energy E of the constructed medium as a function of stages n. The target autocovariance 
function is the Debye random medium function //) given by Eq. (j22p . with a = 30 (pixels) and vol- 
ume fractions (pi = (f)2 = 0.5. The linear size of the system N = 200 (pixels). Surface optimization 
is applied at n = 100, when large clusters have been formed in the constructed medium. 

Numerical experiments show that applying surface optimization properly will further 
decrease the final energy of the (re) constructed medium by a factor of 10~^ (see Fig. [3]), if the 
same coohng schedule is used. Thus, surface optimization enables one to efficiently produce 
more accurate structures associated with the imposed correlation functions. Moreover, by 
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setting a lower energy threshold, one can address the non-uniqueness issue quantitatively. 
Progress on this topic will be reported in our future publications. 

C. Choice of Lattices 

A 2D digitized heterogeneous material (medium) can be represented as a 2D array and 
the real morphology of the material also depends on the lattice on which the system of pixels 
is built . Here, we mainly focus on two commonly used lattices in the literatures, namely, 
the square lattice and the triangular lattice (see Fig. H]). 
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(a) 



(b) 



FIG. 4: (a) The pixels (squares) and intrinsic directions of a square lattice, (b) The pixels 
(hexagons) and intrinsic directions of a triangular lattice. 



Implementation of the algorithms discussed above on a square lattice is easier than that 
on a triangular lattice because the former has orthogonal lattice vectors. However, the pixels 
of a square lattice (squares) only have 4-fold symmetry while those of a triangular lattice 
(hexagons) have 6-fold symmetry, as shown in Fig. HI We find that for the media with 
long-range correlations or large-scale structures, using a square lattice usually introduces 
undesirable anisotropy in the (re) constructed systems, as shown in Fig. [Sj while for the 
media without long-range order or large-scale structures, both lattices work well. 

It is worth pointing out that the triangular lattice is superior to the square lattice in the 
(re) construction of anisotropic materials. In that case, one cannot count on lattice-point 
algorithm, which only uses radial averaged structural information. Several optimization 
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(a) (b) 

FIG. 5: Constructions of media with large-scale structures, (a) Medium generated on a square 
lattice with a square unit cell, (b) Medium generated on a triangular lattice with a rhombus unit 
cell. 

directions need to be specified and processed separately. The triangular lattice has six 
intrinsic directions due to the higher symmetry of its pixel shape while the square lattice 
only has four (see Fig. H]). Methods using more optimization directions have been proposed 



by Torquato and coworkers 



12l |. i.e., along the 45- and 135-degree directions in the square 



lattice; but one needs to pay the cost of complexity of implementation in these cases. 



III. SIGNIFICANCE OF THE ENERGY THRESHOLD 

An important issue that has not been emphasized in our previous work on the 
(re) construction algorithms is the choice of the energy threshold Eth, i.e., the error tol- 
erance of discrepancies between the statistical properties of the generated structure and the 
imposed ones. Recall that in our case, the energy is defined as the sum of squared differ- 
ences between target 5'2(r) and 5'2(r) of the constructed medium [see Eq. ([1])] and the energy 
threshold Eth is a prescribed value of E such that when E < Eth, the (re) construction is 
terminated. Eth indicates how accurately the medium is (re) constructed. A smaller Eth 
means the two-point correlation function of the generated medium matches the imposed one 
better. 

In the following, we will estimate the change of the energy of a digitized medium caused 
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by a perturbation of its original structure (i.e., displacing a randomly selected black pixel). 
This change of energy can be considered as the energy difference between the target and the 
constructed medium. Thus, by imposing a particular value of the energy difference (i.e., Eth), 
we can in turn estimate the number of perturbed black pixels. In this way, we quantitatively 
relate the energy threshold to the question of how well the medium is (re) constructed. In 
general, the energy difference depends on how 5*2 is sampled, the linear size of the system 
and the number of black pixels Nb (volume fraction of black phase ^i), which will be 
discussed accordingly. 



A. Etfi for the Lattice-Point Algorithm 

First, we consider the energy threshold of the lattice-point algorithm. The two-point 
correlation function of a digitized medium is computed from Eq. ([3]). If the generated 
structure perfectly matches the target structure, the energy given by Eq. ([T]) should be 
exactly 0. Suppose the perfect (re) constructed structure is perturbed by moving one of its 
black pixels to an unoccupied lattice site; then 5*2 (r) is different from 5*2 (r) and E becomes 
a small positive number. Note that we assume the perturbed structure does not have the 
same ^2 as the original structure, i.e., there is no structural degeneracy. Let Emm denote 
the smallest positive value of E. From Eq. ([T]), we see that the difference between 5*2 and 
5*2 of every bin contributes to E. From Eq. ([3]) and the fact that the last bin (farthest away 
from the reference center) contains the largest number of pairs of lattice sites, we have 



E„ 
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where [N/2] is the integer part of N/2 and uj^j is the number of elements of the integer set 
Qn, which is defined by 
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n\ 2^ 



(7) 



In other words, Eq. ([6]) estimates Emin using the change of E from caused by removing 
(adding) one pair of black pixels from (to) the last bin (the pixels bounded in pairs due to 
the symmetry of the distance matrix). Also note that the pair of pixels are originally in (or 
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moved to) a "bin" for the pair distance larger than [A^/2], which we do not take into account 
for computation of 5*2. 

For a particular value of Et^, the maximum number of misplaced pairs of black pixels can 
be estimated by 

N^ = p^ = \uj%EtH. (8) 

The ratio of the number of misplaced pairs of black pixels over the total number of pairs of 
black pixels is given by 

Instead of specifying Eth, one can also specify the ratio 7 = 7^; and the threshold can be 
computed from Eq. (Q: 

Eth = -^is(jylN\ (10) 

For example, consider a system composed of 200 x 200 pixels (A^ = 200) and (pi = 0.5. 
The total number of lattice sites is Ns = 40, 000 and the total number of black pixels 
is Nb = 20,000. ti^Ar can be obtained numerically, which is un ~ 10^. From Eq. (EI), 
we have Emm ~ 10~^^. Suppose we choose the threshold Eth = 10~^, from Eq. ([8]), we 
have the number of misplaced pairs of black pixels ~ 10^, which seems to be a large 
number. However, when considering the total number of pairs in the system, we have 
Ntot = = 4 X 10^; from Eq. the ratio is given by 

^ Ntot 4x108 ^ ^ 

which means only one out of a half million pairs is put in the wrong bin. Thus, the medium 
is (re) constructed to a very high accuracy. In our simulations, we choose the threshold 
Eth = 10-^ 

B. Eth for the Orthogonal S'2-Samphng Algorithm 

For the orthogonal S'2-sampling algorithm, 6*2 (r) is sampled line (column) by line (column) 
by moving a line (column) segment of length r one pixel distance each time and counting the 
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times that both ends of the segment are black pixels. This number is then divided by the 
total number of times that one moves the line (column) segment (total number of pixels on 
the line (column)) to obtain 5*2 of that line (column). Finally, the 5*2 sampled from different 
lines and columns of the digitized medium are averaged to compute the 5*2 of the whole 
medium. Similarly, consider the perfect (re) constructed structure is perturbed by moving 
one of its black pixels to an unoccupied lattice site and there is no structural degeneracy, 
is given by 

_ 1 

Emin = (12) 

where N is the linear size of the system. For a particular Eth, the maximum number of 
misplaced black pixels is given by 

N^ = ^ = N'E,,. (13) 

Thus, the ratio of misplaced black pixels over the total number of black pixels can be 
obtained by 

7 = §^ = ^N^E,,. (14) 
Np 01 

If stead, 7 = 7s is specified, the required energy threshold Eth can be obtained from Eq. f|T4l) : 

Consider the same system used in the previous section (A^ = 200, 0i = 0.5). If we require 
only 10 black pixels are misplaced, thus, the ratio 7^ = 10/A^^ = 2.5 x 10~^. From Eq. (fT5l) . 
we have 

2.5 X 10"^ X 0.5 _ „ 

Eth = ^ 10"^ 16) 

4 X 10^ ^ ^ 

In other words, if we choose the threshold Eth = 10~^, only 10 black pixels in the system 
are misplaced. The medium is (re) constructed to a high accuracy. 

It is worth noting that in the above discussions, we assume that the perturbed structure 
does not have the same 5*2 of the original structure, i.e., there is no structural degeneracy. In 
general, however, structural degeneracy does exist (i.e., media with the same 5*2 but different 
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5*3, 5*4, ...)• Results concerning these non-uniqueness issues will be reported in our future 
publications. 



IV. APPLICATIONS OF THE THEORETICAL FORMALISM 

In this section, we illustrate the practical utility of our theoretical formalism in both two 
and three dimensions. In particular, we will consider two kinds of applications. Firstly, 
given a set of basis functions (may not be complete), one can express the scaled autocovari- 
ance functions (two-point correlation functions) of a statistically homogeneous and isotropic 
medium in terms of the specified basis functions with certain accuracy. Realizations of the 
materials can be generated using proper (re) construction procedures and subsequent anal- 
ysis can be performed on the images to obtain effective macroscopic properties of interest; 



see, e.g., Ref. [lOj. Secondly, the map p (see the discussion below) enables one to construct 
candidates of realizable two-point correlation functions with properties of interest, which in 
turn enables one to design and investigate materials with desired structural characteristics. 

The (re) construction of realizations of 3D medium from the information obtained from a 
2D micrograph or image is of great value in practice . Therefore, we will apply the lattice- 
point algorithm to generate three-dimensional digitized realizations of a Fontainebleau sand- 
stone and a boron carbide/aluminum composite from two-dimensional tomographic images 
of their slices through the materials. In a successful reconstruction, other deemed crucial 
structural characteristics besides 5*2 obtained from the reconstructed medium should also 
agree closely with that of the target medium. Consequently, in order to judge quantita- 
tively how well the reconstructions are, we will measure and compare another important 
morphological descriptor, i.e., the two-point cluster function C2(xi,X2), defined to be the 
probability of finding two points at xi and X2 in the same cluster of the phase of interest 4|. 
For statistically homogeneous and isotropic media, C2 only depends on the relative scaler 
distances between the points, i.e., C2(xi,X2) = C2(|xi — X2I) = C2(r). Note that C2 contains 
nontrivial topological "connectedness" information. When large clusters are present in the 
medium, C2 becomes a long-ranged function and its integral will diverge if the phase of 
interest percolates. The measurement of C2 for a 3D material sample cannot be made from 
a 2D cross-section of the material, since it is an intrinsically 3D microstructural function 
To sample C2 from a digitized medium, we associate each pixel with a cluster-index 
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that indicates to which cluster the pixel belongs, and only bin the distances of pixel-pairs 
in the same cluster. Then the number of pair distances in each bin is normalized by the 
total number of distances between the lattice sites within that bin, which is similar to the 
procedure of sampling 5*2 discussed in Sec. II. 

We will also study the reconstruction of a binary laser-speckle pattern in two dimensions 
and show that the algorithm cannot reproduce the pattern accurately. Moreover, we in- 
troduce and discuss a classification of heterogeneous materials into multi-scale media and 
single-scale m,edia. A multi-scale medium (MSM) is the one in which there are multiple 
characteristic length scales associated with different structural elements. Examples of MSM 
include fractal patterns and hierarchical laminate composites 16|. A single-scale medium 
(SSM) is the one composed of structural elements associated with only one characteristic 
length scale, such as a Fontainebleau sandstone and a boron carbide/aluminum composite. 
We conclude that in general reconstructions using 5*2 only work well for heterogeneous ma- 
terials with single-scale structures. However, two-point information via 5*2 is not sufficient 
to capture the key structural features of multi-scale media. Before presenting the aforemen- 
tioned applications, we will first consider a general form of p, which includes all possible 
convex combinations of the basis functions. 



A. A General Form of p 

In paper I, we introduced the idea of expressing the two-point correlation function of a 
statistically homogeneous and isotropic medium through a selected set of bases of the two- 
point correlation function space. In practice, it is convenient to use the scaled autocovariance 
functions that are equivalent to the two-point correlation functions, which are defined as 
follows: 

/(.) . (17) 

0102 

Suppose is a set of bases of scaled autocovariance functions and p is an map on 

{/i(r)}™ composed of convex combinations and products of {i = 1, ...,m), thus 

fir) = p[{/.(r)}r=i] = p[/i(r), /2(r), /„(r)], (18) 
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is also a realizable scaled autocovariance function. Note that the choice of basis functions is 
not unique. 

In general, one can consider that takes the form 

Pmr)}T=i] = E«^/^W + E A,/.(r)/,(r) + • • • , (19) 
where the coefficients satisfy the condition 

E«^ + EA, + --- = i. (20) 

Then one can use standard regression methods to obtain the set of coefficients such that 
/(r) = p[{/j('")}i^i] is the best approximation of the target function pointwisely. 

The basis functions we consider here include Debye random medium function /^(r), the 
family of polynomial functions fp'\r), damped oscillating function fo{^), overlapping-sphere 
function fs{f) and symmetric-cell material function fc{'^)i which are discussed in paper I. 



B. Modeling Real Materials 

1. Fontainebleau Sandstone 



First, we investigate structural properties of a Fontainebleau sandstone from a two- 
dimensional tomographic image of a slice through the material sample. Sandstone is an 
important porous medium in g eo-p hysical and petroleum applications and has been the fo- 



cus of many studies 



13, 



1^ 



191]. A microstructural image of a slice of a Fontainebleau 



sandstone is shown in Fig. El^a), in which the black solid phases (phase 1) and the 

white areas are void phases (phase 2). The two-point correlation function of the void phase 
S2^{r) is shown in Fig. [6](b). 

The scaled autocovariance function of the void phase in Fontainebleau sandstone f^^{r) 
can be approximated by the convex combination of /^(r) and fo{f) as follows: 



f^{r) = aMr) + a2fo{r), 



(21) 



where ai = 0.77, a2 = 0.23 are the combination coefficients and 
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(a) 



40 
r (pixels) 



(b) 

FIG. 6: (a) A microstructrual image of a slice of a Fontainbleau sandstone The black phase 
is the solid phase with = 0.825, and the white phase is the void phase with (j)2 = 0.175. (b) The 
two-point correlation function 82'^ (r) of the void (white) phase. 



foir) = exp(-r/a), 



(22) 



fo{r) = exp(— r/6) cos(gr + ip), (23) 

where a = 3 [pixel), b = 6.5 (pixel) are the effective correlation length; q = 0.2 {pixel~^) is 
the oscillating frequency and ip = is the phase angle. The two-point correlation function 
S[^{r) is approximated by 
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(a) (b) 

FIG. 7: (color online), (a) The reconstructed 3D realization of the Fontainbleau sandstone from 
(r) of the void phase. The void phase is shown in yellow and the solid phase is transparent for 
easy visualization, (b) A 2D slice of the constructed 3D realization. The solid phase is shown in 
black and the void phase is shown in white. The linear size of the system N = 160 (pixels). 

S^ir) = f'{r)M2 + 0^ + 5^2(r), (24) 

where (pi = 0.825, 02 = 0.175 are volume fractions of the solid and void phase, respectively; 

(r) is the discrepancy between the sampled two-point correlation function and the basis- 
function approximation. The average of the absolute values of discrepancies AS'2, which 
indicates how well the sampled two-point correlation function is approximated by the convex 
combination, is defined as: 

^•^2 = i^El'^'^2(r)| ~ 5.2 X 10-^ (25) 

where = 80 is the sample-length. Note that although the general form of p given by 
Eq. ( JTOl) would work well if a complete set of basis functions is given, in the present case, 
our practical approach is enough. 

The 3D reconstruction of the Fontainebleau sandstone from S'^'^(r) obtained from the 
digitized image of a 2D slice [Fig. [6](a)] is shown in Fig. [3 Visually, the reconstruction pro- 
vides a good rendition of the true sandstone microstructure, as can be seen by comparing the 
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FIG. 8: The two-point cluster functions (r) for the void (white) phase of the target and recon- 
structed 2D shces of the Fontainebleau sandstone. 

2D images. As pointed out earlier, to ascertain whether the reconstruction is quantitatively 
successful, we also measure C2 of the target and generated media. Since it is an intrinsically 
3D microstructural function, we only compute and compare the two-point cluster functions 
of the 2D slices. The measured C2{r) for the void phase of the target and the reconstructed 
slices of the Fontainebleau sandstone are shown in Fig. [HI The figure reveals that although 
C2{r) of the generated medium is slightly below that of the target medium, the discrep- 
ancies are acceptable (e.g., the largest discrepancy I^Cal ~ 5 x 10~^). Thus, we consider 
the reconstruction is successful. Also note that the close agreement of the two-point cluster 
function indicates that C2{r) is largely determined by (r) of the medium and suggests 
that C2{r) might be expressed as a functional of S'^'^(r). 

2. Boron Carbide/ Aluminum Composite 

A 2D digitized image of a boron carbide/ aluminum (-B4C / Al) inter-penetrating composite 
and the sampled two-point correlation function of the aluminum phase (white phase) 5'^'(r) 
{20! are shown in Fig. [91 As we can see from Fig. Mjo), 5'^'(r) is essentially an exponentially 
decreasing function without any significant short-range correlation. Thus, S'^'(r) can be 
approximated by 

S^\r) = a[M2fD{r) + a^0i02/o(r) + 0^ + 55^(r), (26) 
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(a) 
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(b) 



FIG. 9: (a) A digitized image of a boron carbide/aluminum composite 20(]. The black phase 
is boron carbide with = 0.647, and the white phase is aluminum with (p2 = 0.353. (b) The 
two-point correlation S'2^'(r) of the aluminum phase. 



where a'^ = 0.81, ctg = 0.19 are combination coefficients; = 0.647 and 02 = 0.353 
are volume fractions of the boron carbide (black) phase and the aluminum (white) phase, 
respectively, foir) and fo{r) are given by Eq. (1221) and Eq. fl2^ . respectively; the parameters 
used are a = 3 (pixel), 6 = 10 (pixel) and q = 0.22 {pixel~^), = 0. The average of the 
absolute values of discrepancies AS'2 is given by 



^^2 = ^EI'^^2WI^ 6.4x10-1 



(27) 
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(a) (b) 

FIG. 10: (color online), (a) The reconstructed 3D realization of the boron carbide/aluminum com- 
posite from 5*2^' (r). The aluminum phase is shown in yellow and the ceramic phase is transparent 
for easy visualization, (b) A 2D slice of the constructed 3D realization. The ceramic phase is 
shown in black and the aluminum phase is shown in white. The linear size of the system = 160 
(pixels) . 




FIG. 11: The two-point cluster functions C2 '(^") for the aluminum (white) phase of the target and 
reconstructed 2D slices of the boron carbide/aluminum composite. 
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where Nl = 80 is the sample-length. 

The 3D reconstruction of the boron carbide/aluminum composite from S'^'(r) obtained 
from the digitized image of a 2D slice [Fig. M^a)] is shown in Fig. [TUl As in the previ- 
ous section, to verify that the reconstruction is quantitatively successful, we compute and 
compare the two-point cluster function C^'(r) for the aluminum phase of the target and 
reconstructed 2D slices of the composite. As shown in Fig. [TTl the reconstruction also 
slightly underestimates C^'(r) with acceptable discrepancies (e.g., the largest discrepancy 
|(5C^'| ~ 8 X 10^'^). The medium can be considered as successfully reconstructed. Also note 
that C^'(r) is long-ranged, indicating large clusters of the aluminum phase present in the 
media. 

3. Binary Laser-Speckle Pattern 

Although in the above two examples our theoretical formalism works well, there are 
situations where the microstructural information contained in 5*2 (r) alone is not sufficient 
to accurately reconstruct a heterogeneous material. One such example is the multi-scale 
structure of a binary laser-speckle pattern (Fig. [T2l) . The figure reveals that there are three 
structural elements: "particles" , "stripes" and a background "noise" (individual black pixels 
dispersed throughout the white phase). Thus, there are three characteristic length scales in 
the medium associated with these structural elements. 

The reconstruction of the speckle pattern is shown in Fig. [131 Comparing Fig. [13] with 
Fig. [121 we can see that instead of reproducing all the structural elements in the target 
medium, the (re) construction program seems to mix them up to generate a single-scale 
structure that has the same (or to a very high accuracy) two-point correlation function as 
the target medium. We note that this is a numerical example of structural degeneracy of 
S2{r). 

The laser-speckle pattern we considered is also an example of the multi-scale media 
(MSM). The separation of length scales in MSM results in the inefficiency of the bulk- 
based structure characteristics (e.g., ra-point correlation functions), since usually they can 
only pick up structural information associated with the largest length scale. For example, 
consider the dislocations in a crystalline solid as the phase of interest; it is clear that the 
two-point correlation function of the "dislocation" phase is identically zero since the "dislo- 
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FIG. 12: (a) A digitized image of a binary laser-speckle pattern |2l|. The volume fraction of the 
black phase is = 0.639, and the volume fraction of the white phase is (j)2 = 0.361. (b) The 
two-point correlation 5f^(r) of the black phase. 

cation" phase has no measure compared with the bulk "crystal" phase (i.e., the volume of 
dislocations is zero compared with that of the bulk crystal). In this extreme example, the 
ratio 7 of the two characteristic length scales associated with the "dislocation" phase Ld and 
the "crystal" phase is zero, i.e., 7 = L^/Lc = 0. In digitized media, the ratio of length 
scales is always positive due to the discrete nature of the system; however when the ratio 
is significantly different from unity, the two-point correlation function alone is not able to 
capture the key structural features. On the other hand, extensive experience with successful 
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FIG. 13: The reconstruction of the binary laser-speckle pattern from (r). The hnear size of 
the system = 150 (pixels). 



reconstructions of single-scale media (SSM) from 5*2 shows that the two-point correlation 
functions are indeed sufficient to determine the structures of SSM to a high accuracy. Thus, 
we conclude that in general reconstructions using 5*2 work well for heterogeneous materials 
with single-scale structures, while the microstructural information contained in S2 is not 
sufficient to accurately model multi-scale media. Note that more morphological informa- 
tion (e.g., the lineal-path function 221] and the two-point cluster function, etc.) can be 
used to model MSM more accurately. However, even the cluster-type functions containing 
connectedness information of the media such as C2 can not completely characterize MSM 
statistically. 



C. Constructing Materials with Desired Structural Characteristics 

From Eq. (fT8|) . one can construct candidates of realizable two-point correlation functions 
using the basis functions. Given a set of basis functions with diverse and interesting proper- 
ties, one would be able to construct a two-point correlation function that exhibits all the use- 
ful properties of the basis functions to some extend and generate an "optimal" structure that 
realizes all the desired structural features. Thus, the theoretical formalism enables one to 
design materials with structural properties of interest. Given an accurate structure-property 
relation, one could even design materials with physical properties of interest by manipulat- 
ing their two-point correlation functions. For example, Adams, Gao and Kalidindi recently 
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developed a methodology to obtain finite approximations of the second-order properties clo- 
sure in single phase polycrystalline materials, from which the second-order microstructure 
design can proceed [23!. 

As pointed out in paper I, we know very little about the basis function set {/j(r)}™^ at 
this stage. Our choice of is based on the criteria that /i(r)'s should have simple 

analytical forms and they present typical features of certain known two-point correlation 
functions. Important features exhibited by most realizable two-point correlation functions 
are monotonically deceasing or damped oscillating, which corresponds to materials without 
or with significant short-range order, respectively. Another feature of two-point correlation 
functions that may affect the structures of the corresponding materials is the smoothness 
of the function, which is not emphasized in our previous work. In the following, we will 
see through several examples that the properties of basis functions can be observed in the 
generated structures; also the media with desired structural properties can be obtained 
by manipulating the combination coefficients and the parameters (e.g., effective correlation 
length) in the basis functions. 

In the paper I, we already investigated hypothetical correlation functions combining the 
exponentially decreasing and damped-oscillating features. In this section, we provide an 
example of non-smooth correlation functions (with discontinuous derivatives), i.e., the poly- 
nomial function of order two fp\r), which is defined as 



'1 — r/c)^ < r < c, 

(28) 

r > c, 

where the parameter c is the effective correlation length. Other correlation functions hav- 
ing this non-smooth feature include the functions of known constructions (e.g., /s'(r), fc{r), 
etc.) and other polynomial functions 2I. A typical medium associated with this type of func- 
tions is composed of dispersions of fully penetrable "particles" . The monotonically decaying 
part of the functions determines the size and shape of the "particles" and the penetrability 
is consistent with the flat "tail" of the functions, which implies no spatial correlation be- 
tween the "particles". When these functions are use as dominant basis functions, one can 
expect the generated media also contain well-deflned overlapping "particles"; thus, surface 
optimization could be applied in the (re) constructions. However, when functions like Debye 
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random medium function are dominant in the convex combination, this algorithm modi- 
fication should be used with care because the media could contain "clusters" of all sizes 
and shapes, which would significantly reduce the efficiency of the modified algorithm using 
surface optimization. 

Besides fp \r), we also use Debye random medium function fDij") and damped-oscillating 
function /o('") as the basis functions. Consider the simplest form of p for these three basis 
functions, i.e., 

/(r) = «i^(r) + a2/o(r) + «3#(r), (29) 

where a, {i = 1,2,3) satisfies < < 1 and J^i^a = 1; and fnir), fo{r) and fp {r) 
are given by Eq. (!22|) . Eq. (!23|) and Eq. (128!) . respectively. The characteristic length scales 
of the generated structures are determined by the parameters in the basis functions; and 
the ratios of the characteristic lengths are chosen to be close to unity, i.e., the hypothetical 
media belong to SSM. 
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FIG. 14: The constructed autocovariance function /(r) given by Eq. (j29p . The combination 
coefficients are (ai, 02, 03) = (0.3, 0.2, 0.5); the parameters in the basis functions are a = 8 {pixel); 
6 = 5 {pixel), q = 1 {pixel)~^, -0 = 0; c = 5 {pixel). 

Suppose we choose the combination coefficients {ai,a2,a3) = (0.3,0.2,0.5); and choose 
the following values for the parameters in the basis functions: a = 8 {pixel); 6 = 5 {pixel), 
q = 1 {pixel)~^, ip = 0; c = 5 {pixel). The combined autocovariance function f{r) is 
shown in Fig. [14] and the constructed structure with volume fraction of black phase 0i = 0.5 
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(a) (b) 

FIG. 15: (color online), (a) The constructed 3D structure associated with /(r) shown in Fig. [TJ] 
with volume fractions (pi = 02 = 0.5. The "black phase" is shown in yellow and the "white phase" 
is transparent for easy visualization, (b) A 2D slice of the constructed 3D realization. The linear 
size of the system = 150 (pixels). 

is shown in Fig. [T3 From the figures, we can see the effects of each basis function on the 
generated structure. The short-range correlations in the structure are determined by fp (r), 
which allows spatially uncorrelated "particles" with diameter c to form. These "particles" 
would form clusters of different sizes as the volume fraction of the black phase increases. 
The middle-range correlations in the structure are dominated by the oscillation part in 
/(r), which is the contribution of foi^)- As we have shown in paper I, the parameter q 
is manifested as a characteristic repulsion among different structure elements (appearing in 
different forms at different volume fraction) with diameter of order q; while b controls the 
overall exponential damping, and thus, the effective range of the repulsion. The large-scale 
correlation is dominated by the long "tail" of f 0(1^)7 which allows clusters of all shapes 
and sizes to form. In Fig. [15], we can clearly identify the structure elements associated 
with each basis function. For example, the stripe-like structures are associated with the 
oscillating feature of f{r), the width of which is approximately q. Several "particles" with 
diameter a can be found dispersed in the white phase, which correspond to the contribution 
of fp (r). Note that most "particles" form clusters together with the "stripes" at this 
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relatively high volume fraction 0i = 0.5. Finally, the spatial distribution of large-scale 
clusters is determined by the "tail" of /z)(r). 
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FIG. 16: The constructed autocovariance function /(r) given by Eq. (I29p . The combination 
coefficients are (ai, 02, as) = (0.3, 0.2, 0.5); the parameters in the basis functions are a = 8 (pixel); 
6 = 10 (pixel), q = 0.5 (pixel)~^, 'ip = 0; c = 15 (pixel). 

Suppose now we would like to generate a similar structure as the one in Fig. [TS] but with 
larger "particles" and clusters. To "grow" the "particles", we need to increase the param- 
eter c in fp\r), which controls the effective diameter of the "particles". To form larger 
clusters, we need to reduce the "repulsion" between the structural elements; or equivalently, 
to suppress the oscillation introduced by /o(^)- This can be done by either increasing the 
parameter b in fo{r) or decreasing the coefficient a2- We first adopt the former method. 
Thus, to construct a new structure with the required properties, we use the following pa- 
rameters: a = 8 (pixel); 6 = 10 (pixel), q = 0.5 (pixel)^^, ■?/' = 0; c = 15 (pixel); and use 
the same (i = 1,2,3) as the last example. The combined autocovariance function f(r) is 
shown in Fig. [16] and the constructed structure with volume fraction of black phase 0i = 0.5 
is shown in Fig. [171 From the figures, we can see the generated medium indeed contains 
larger "particles" and clusters; besides, there are still stripe-like structures associated with 
the repulsion part of f(r). 

Now we would like to further increase the size of clusters in the constructed medium while 
keeping all the parameters in the basis functions unchanged. This time we need to adjust 
the combination coefficients to obtain the required structure. To generate larger clusters. 
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FIG. 17: (color online), (a) The constructed 3D structure associated with /(r) shown in Fig. [16] 
with volume fractions (pi = 02 = 0.5. The "black phase" is shown in yellow and the "white phase" 
is transparent for easy visualization, (b) A 2D slice of the constructed 3D realization. The linear 
size of the system = 150 (pixels). 

the repulsion effect need to be further suppressed; thus, we need to reduce ^2. Since the 
sum of the coefficients must equal 1, we also need to increase the other two (z = 1,3). 

Note that further increasing 0:3 makes fp {r) dominant in the combination; however, as 

(2) 

we discussed before, fp {r) corresponds to spatial uncorrelated "particles" and the size 
of clusters formed by these "particles" is determined by the volume fraction the "particle" 
phase, which cannot be changed here. So we choose to increase ai of /^(r) based on the fact 
that the large correlation length and long "tail" of fnif') significantly stimulate the forming 
of large clusters. Thus, in the construction we use the following combination coefficients: 
(«!, a2, 0:3) = (0.45, 0.05, 0.5); and use the same parameters in the basis functions as the last 
example. The combined autocovariance function /(r) is shown in Fig.[T8]and the constructed 
structure with volume fraction of black phase 0i = 0.5 is shown in Fig. [191 

From the above examples, we see that even a simple convex combination of the basis func- 
tions (the simplest form of map p) can provide variety of candidates of realizable two-point 
correlation functions. We can also obtain desired structures by manipulating the combina- 
tion coefficients and the parameters in the basis functions. Thus, we argue that in general 
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FIG. 18: The constructed autocovariance function /(r) given by Eq. (j29p . The combination coef- 
ficients are (01,02,03) = (0.45,0.05,0.5); the parameters in the basis functions are a = 8 (pixel); 
6 = 10 (pixel), q = 0.5 (pixel)^^, ip = 0; c = 15 (pixel). 

given a complete set of basis functions, Eq. (fT9l) enables one to construct candidates of re- 
alizable functions with desired properties and to design materials with structural properties 
of interest. 



V. CONCLUSIONS 



In this paper, we described in great detail the lattice-point algorithm, which has been 
shown to be both efficient and isotropy-preserving for (re) constructing statistically homo- 
geneous and isotropic media. The digitized medium (pixel system) can be considered as a 
lattice-gas system, in which the black pixels behave like nonoverlapping "particles" moving 
from one lattice site to another. The two-point correlation function of the medium is sampled 
by binning all distances between black pixels and dividing the number in each bin by the total 
number of distances between all lattice sites in that bin. The energy threshold indicating to 
what accuracy the medium is (re) constructed has been discussed in detail for the first time. 
This quantity is directly related to the issue of non-uniqueness of (re) constructions. We 
have also described an algorithm modification using surface optimization to further speed 
up the (re) construction process and discussed its implementation based on our lattice-gas 
version of the digitized media. Numerical experiments have shown that by applying the 
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(a) (b) 

FIG. 19: (color online), (a) The constructed 3D structure associated with /(r) shown in Fig. [18] 
with volume fractions (pi = 02 = 0.5. The "black phase" is shown in yellow and the "white phase" 
is transparent for easy visualization, (b) A 2D slice of the constructed 3D realization. The linear 
size of the system = 150 (pixels). 

surface optimization algorithm properly (i.e., in the (re) constructions of media composed 
well-defined "particles"), the final energy can be reduced by a factor of 10^^ if the same 
cooling schedule is used. The choice of different pixel-lattices has also been discussed. 

We have applied the theoretical formalism proposed in paper 1 to model several exam- 
ples of real materials. In particular, we used the lattice-point algorithm to generate 3D 
realizations of a Fontainebleau sandstone and a boron carbide/aluminum composite from 
2D images of their slices through the materials. The two-point cluster functions in the re- 
constructions in these two different cases matched the corresponding ones in the original 
materials, thus demonstrating the efficiency and accuracy of the (re) construction algorithm 
in these cases. We also studied the reconstruction of a binary laser-speckle pattern in 2D, 
in which the algorithm fails to reproduce the target pattern accurately. We conclude that 
in general reconstructions using 5*2 only work well for heterogeneous materials with single- 
scale structures. However, two-point information via 5*2 is not sufficient to accurately model 
multi-scale media. In paper I, we pointed out that given a complete set of basis functions, 
one can obtain the basis-function approximations of the two-point correlation functions for 
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the materials of interest to any accuracy. Here we only used a simple form of Eq. (fT9l) due 
to our limited knowledge of the basis functions set. We also constructed realizations of ma- 
terials with desired structural characteristics by manipulating the combination coefficients 
and the parameters in the basis functions. Thus, Eq. f[T^ enables one to construct candi- 
dates of realizable functions with desired properties and to design materials with structural 
properties of interest. 

We are now developing efficient (re) construction procedures that take into account ad- 
ditional microstructural information that characterize the media. For example, using our 
lattice-gas version of the digitized media, the two-point cluster functions C2 can be directly 
incorporated into the Yeong-Torquato scheme 5j], which enables one to obtain more accu- 
rate (re) constructions of the target media. Also, a quantitative structure-property relation 
is necessary for the further applications of the theoretical formalism. Finite element anal- 
ysis on the generated media will be performed to establish the relation. Such work will be 
reported in our future publications. 

Acknowledgments 

Acknowledgment is made to the Donors of the American Chemical Society Petroleum 
Research Fund for support of this research. 



[1] S. Torquato, Random Heterogeneous Materials: Micro structure and Macroscopic Properties 

(Springer- Verlag, New York, 2002). 
[2] Y. Jiao, F. H. Stillinger, and S. Torquato, Phys. Rev. E 76, 031110 (2007). 
[3] C. L. Y. Yeong and S. Torquato, Phys. Rev. E 58, 224 (1998). 
[4] S. Torquato, J. D. Deasley, and Y. C. Chiew, J. Chem. Phys. 88, 6540 (1988). 
[5] C. L. Y. Yeong and S. Torquato, Phys. Rev. E 57, 495 (1998). 

[6] K. Wu, M. I. J. Dijke, G. D. Couples, Z. Jiang, J. Ma, K. S. Sorbie, J. Crawford, I. Young, 

and X. Zhang, Trans. Porous Media 65, 443 (2006). 
[7] M. A. Ansari and F. Stepanek, AlChE Journal 52, 3762 (2006). 
[8] R. Hilfer and C. Manwart, Phys. Rev. E 64, 021304 (2001). 



34 



[9] D. Basanta, M. A. Miodownik, E. A. Holm, and P. J. Bentley, Metall. Mater. Trans. A 36, 
1643 (2005). 



[lo: 

[11 
[12: 
[13 
[14 
[15: 
[16 

[17; 

[18 

[19 
[20 

[21 
[22; 
[23 



H. Kumar, C. L. Briant, and W. A. Curtin, Mech. Mater. 38, 818 (2006). 

D. Cule and S. Torquato, J. Appl. Phys. 86, 3428 (1999). 

N. Sheehan and S. Torquato, J. Appl. Phys. 89, 53 (2001). 

M. G. Rozman and M. Utz, Phys. Rev. Lett. 89, 135501 (2002). 

R. Hosemann and S. N. Bagchi, Acta Cryst. 7, 237 (1954). 

C. Chubb and J. I. Yellott, Vision Research 40, 485 (2002). 

J. Quintanilla and S. Torquato, Phys. Rev. E 53, 4368 (1996). 

D. A. Coker, S. Torquato, and J. Dunsmuir, J. Geophys. Res. 101, 17497 (1996). 

L. M. Schwartz, F. Auzerais, J. Dunsmuir, N. Martys, D. P. Bentz, and S. Torquato, Physica 
A 207, 28 (1994). 

C. Manwart, S. Torquato, and H. R., Phys. Rev. E 62, 893 (2000). 

S. Torquato, C. L. Y. Yeong, M. D. Rintoul, D. L. Milius, and I. A. Aksay, J. Am. Ceram. 
Soc. 82, 1263 (1999). 



The source of the original gray-scale image: http://www.people.vcu.edu/~ecarpenter2 



B. Lu and S. Torquato, Phys. Rev. A 45, 922 (1992). 

B. Adams, X. Gao, and S. Kalidindi, Acta Mater. 53, 3563 (2005). 



35 



